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Abstract. We analyze a single-level quantum system placed between metallic 
leads and strongly coupled to a localized vibrational mode, which models a single- 
molecule junction or an STM setup. We consider a polaron model describing the 
interaction between electronic and vibronic degrees of freedom and develop and 
examine different truncation schemes in the equation-of-motion method within the 
framework of non-equilibrium Green functions. We show that upon applying gate 
or bias voltage, it is possible to observe chargc-bistability and hysteretic behavior 
which can be the basis of a charge-memory element. We further perform a 
systematic analysis of the bistability behaviour of the system for different internal 
parameters such as the electron-vibron and the lead-molecule coupling strength. 
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1. Introduction 

Within the field of single- molecule electronics [2 El E] , beside experimental progress 
with regard to vibrational properties and their signatures in transport J5j El [9] , 
related phenomena such as switching, memory effects and hysteretic behavior in 
molecular junctions have gained increasing importance and attention. Random and 
controlled switching of single molecules [TDl [TTJ [12] , as well as conformational memory 
effects [T31 HH [H] have been recently explored. Related to these effects, there is the 
so-called charge-memory effect, that is basically a hysteretic behaviour of the charge- 
voltage, respectively, current-voltage characteristics arising from the interplay between 
the polaron shift and Franck-Condon blockade [16j in the presence of electron- vibron 
interaction. Several works have recently addressed this interesting feature of molecular 
junctions, both experimentally [12 EH] and theoretically [19 l [20 l [2l l [22 [ [23] . 

The charge-memory effect can be explained in the framework of a simple single- 
level polaron model [20] E31 EH ES] , where the electronic state is coupled to a vibronic 
mode with frequency u>o, see the sketch in Fig. [TJ If the energy of the unoccupied 
electron level without electron- vibron interaction is eo ; the occupied (charged) state 
of the interacting system will have the energy ei = to — e p , where e p is so-called 
polaron shift (called also recombination energy). Neutral and charged (polaron) states 
correspond to local minima of the potential energy surface and get mctastablc, if the 
electron-vibron interaction is strong enough. Upon applying an external voltage, one 
can change the state of this bistable system, an effect that is accompanied by hysteretic 
charge- voltage and current-voltage curves. In this model it is not necessary to include 
Coulomb interaction explicitly, though one can straightforwardly incorporate charging 
effects. 

It was alternatively suggested in Refs. E2 E2] that quantum switching between 
bistable states rather results in telegraph noise at finite voltage than in a memory 
effect. In a recent paper [23], considering the problem in the weak molecule-to-lead 
coupling limit, i.e. for the level width Y <C luq, e p , we considered the crossover between 
these two pictures, if one takes into account the time scale of the switching process. 
Indeed, the switching time t between the two states of interest should be compared 
with the characteristic time of the external voltage sweeping, r s ~ V (t) / {dV (t) / dt) . 
For r 3> t s , quantum switching can be neglected, and hysteresis can be observed, 
while in the opposite limit, r <C t s , the averaging removes the hysteresis. In Rcf. [23] 
we found that at large enough electron-vibron interaction strcnght \/u>q (where A is 
introduced below in Eq. {!])), the switching time r is exponentially large compared to 
the inverse bare tunneling rate 1/r. 

In this paper we consider the polaron problem in another limiting case, luq <C Y < 
6 p , corresponding to the regime where the Born-Oppcnhcimcr approximation holds 
true. In this situation of intermediate molecule-to-lead coupling we address the case 
of time-independent applied voltages, considering then the stationary problem and 
focusing on the properties of the two states of interest. In this paper we consider the 
parameter ranges for which the fluctuations between the two charge states of interest 
are negligible. It should be noted that the mean-field solutions we consider are meta- 
stable, but they are physical in the case of very long switching time r. Of course, 
the exact equilibrium solution is a superposition of these two states. We assume 
that all possible relaxation processes, with the exception of the switching between 
metastable states, are much faster than r. Within the framework of nonequilibrium 
Green functions [26l [27l [28l [29l [30] we use an equation-of-motion (EOM) approach 
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Figure 1. (Color) Schematic representation of a single-level model system 
interacting with a vibronic mode and coupled to left and right leads. 



which allows for studying the appearance of the charge-memory effect addressed at 
different levels of approximation, starting with the self-consistent Hartree level. The 
work is partly built on and further develops ideas introduced in Refs. [24] ESJ [20] and 
complements the results based on a master-equation approach in Ref . [55] . 

The paper is organized as follows: In Sec. 2 we introduce the model and the 
relevant quantities necessary to describe transport through the molecular junction. 
In Sec. 3 we outline the equation-of-motion method, used for the calculation of 
nonequilibrium Green functions, and discuss different truncation schemes. In the 
final Sec. 4 we present and discuss the results of our approach. 

2. Model and method 

2.1. The single-level electron-vibron Hamiltonian 

We consider a three-terminal device with the central part given by a single level 
interacting with a vibronic mode (sec Fig. [I]). Possible experimental implementations 
include metal-molecule-metal junctions or STM spectroscopy of a single molecule on a 
conducting substrate. To describe such a system we use the standard electron-vibron 
Hamiltonian 

+ A(a f +a)dU + 
+ [( £i > fc + e< ^) c U c a + ( V hk4,k d + h - c -) ■ (!) 

i.k 

Here the terms in the first line describe the central system including one electronic 
state with energy eo, one vibronic state with frequency too and their mutual interaction 
with a coupling strength A. The second line in Eq. (fTJ) contains the Hamiltonian of 
the two leads with independent-particle states and the tunneling between the leads 
and the central region via the couplings V^k- The index i denotes the left and right 
leads, while k labels the electronic states of electrons in the leads. Below it proves 
furthermore convenient to employ the vibronic " position" and " momentum" operators 

x = (V + a ; p = — a . (2) 
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The energy level eo in Eq. (|TJ) can be shifted through the gate voltage Vq. 
We choose as reference energy e = for Vq = and assume a linear capacitive 
coupling. e = aeVc putting a = 1. Note that the presence of a bias voltage 
Vb — <Pl — <Pr can also change the energy of the electronic level via the parameter 

= PR + WVb, where < 77 < 1 describes the symmetry of the voltage drop 
across the junction: rj = corresponds to the completely asymmetric case, while 
rj = 0.5 stands for the symmetric case. Thus in all approximations the bias and gate 
voltages are taken into account twofold, through the potentials of the leads, which 
can be chosen e.g. as <pl = Vb/2, tf R = —Vb/2, and the effective energy of the 
level, which is correspondingly e = eo + e( Po = eVc + e[r\ — 0.5)Vb- From this 
expression for eo it follows that in the case of asymmetric bias-voltage drop across 
the junction, 77 = 0, the energy of the unoccupied electron level will be centered 
around the electrochemical potential of the right lead and moved away from this value 
through the gate voltage. This ingredient will be crucial for the effect addressed 
in this paper because the additional presence of the polaron shift will then fix the 
energy of the occupied (charged) state below the electrochemical potential of the right 
lead resulting in a blocked charged state under appropriate parameter conditions. In 
the case of symmetric voltage drop, r) = 0.5, the energy of the unoccupied electron 
level will be centered around the zero of the energy resulting in a different scenario for 
what concerns the memory effect. In this paper we consider the case of an asymmetric 
junction for which we show that the memory effect occurs for small bias voltage in 
a wide range of the parameters entering the model Hamiltonian (JlJ. The symmetric 
situation can also give rise to a hysteretic behaviour but only at finite bias voltages, 
making this case less interesting for the memory effect addressed here. 



2. 2. Spectral function, average charge and current 

To obtain physical information from the Hamiltonian {!]) we use the Green function 
method within the equation of motion (EOM) formalism. This method is an 
alternative to the Green function techniques earlier applied to the considered problem 
in Refs. [20l [23l EH [25]. We start with the retarded Green function for two generic 
operators a and b which is defined as 

G r (t u t 2 ) = -i& (h -t 2 ) ({a^Xh)}) , (3) 

where (, ) denotes a thermal average and {a, b} — ab + ba is the anti-commutator. For 
the stationary case, the general expression © reduces to an object with only one time 
argument, 

G r (t) = -te(t)({a(t),b}) ■ (4) 
For a single-level system the retarded Green function in Eq. @ is obtained by replacing 
the generic operators a and b by the electronic operators. For the Hamiltonian |T|), 
this reads 

G r (t) = -ie(t){{d(t),J}) . (5) 

From Eq. ([5]) one obtains the spectral function A(e) of the system through the 
expression 

A(e) = -2ImG r (e), (6) 

where G r (e) is the Fourier transform of G r (t). The spectral function is the basic 
ingredient for obtaining the transport properties of the system such as average current 
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and charge on the molecule. The expression for the current through the molecule is 
given by 

/= T?¥^- f~^)[/£(e-e^)-/°(e-e^)]|i, (7) 

where is the equilibrium Fermi function in the z-th lead. The tunneling couplings 
to the right (T R ) and left (T L ) leads are 

Ti(e) = 2^^|^ fe | 2 5(e-e 4 , fc ), (8) 

k 

where the matrix elements Vi t k are assumed to be energy-independent (wide-band 
limit). The full level broadening is given by the sum r = Tl + Tr. Below r# and 
are assumed to be the same. 

The average charge (number of electrons), n = (d'd\ is given by 

+oc de 

A{e)f{e)-, (9) 

-oo 2lr 

where /(e) is the distribution function of electrons inside the molecule. For the 
approximations used in this paper, we employ the same distribution function as in 
the non-interacting case, 

/W-Esfife^lEsfife^, (10) 

1 L + 1 R 

because we are focusing on the case of intermediate moleculc-lcad coupling. Fast 
tunneling into and out of to the molecule makes plausible the assumption that the 
electrons are in strong a non-equilibrium situation and can then be decribed via 
Eq. (fTU]) that is obtained assuming only elastic processes. Moreover, within the first 
approximation that we discuss in the paper, it is possible to obtain the distribution 
function (|10[) analytically by calculating the lesser Green function of the problem or 
applying the Hartree approximation directly to the Hamiltonian as in Ref. |20j . 

For a more complete explanation of the basic formulas introduced here the reader 
is referred to Refs. EE1 ESI [30] . 



3. Equation-of-motion method 

3.1. General formalism 

The equations of motion for non-equilibrium Green functions arc obtained from the 
Heisenberg equation for a Heisenberg operator a(t), 

.da 



a,H 



aH-Ha. (11) 



Here and below all Hamiltonians are assumed to be time-independent because we 
consider the stationary problem, the applied voltages enter the problem only as time- 
independent parameters, changing the position of the molecular energy level and the 
electrochemical potentials in the left and the right lead. With Eq. (|TTj) the time 
derivative of the Green function, Eq. (j4|), reads 

i—^J- = S(t) ({a(t), b}) - z6 (f) ({[a, H] (t), b}) . (12) 
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After performing a Fourier transform of Eq. Q12p we obtain 

(e + irj){{a,b)) = ({a,b}) + {([a,H],b)) , (13) 

where ((, }) indicates the Fourier transform of a given Green function. Equation (|13[) 
is the starting point for the EOM method. By applying successively a time derivative 
there are new high-order Green functions appearing. The idea of the EOM method 
is thereby to truncate this iterative process at some point making the mean-field like 
approximation of the highest-order Green function through lower-order functions, in 
order to obtain a closed set of equations. In our case we start from the Hamiltonian 
|T]) and consider the first-oder equation for the function ((d,a f }} and second-order 
equations for the functions ((xd, d^\) and ((pd, d^)). 

3.2. EOM method for the single-level electron-vibron Hamiltonian 
The EOM method for Hamiltonian |T]) generates the expressions 

(e + ty) ((<*,«**» =l + eo((d,d^))+X((xd,d^)) 

+ £^« c ^' dt )) , (14) 

(e + in- e i>k ) = V itk ((d, dt» . (15) 

The equation for ((ci^jdn) is closed (including only the function ((d, d')))- By 
substituting Eq. (fTT))) into Eq. (fl"4|) and introducing the self-energy E of the leads 
through 

e = v ^ fc|2 , (i6) 

we obtain eventually 

(e + in - e - E) {{d, tfi)) = 1 + \((xd, d f )) . (17) 

The last term, describing the interaction between electron and vibron, has to be 
truncated at this level or found from higher-order equations and then truncated at a 
higher level of approximation. The lead self-energy will be used below in the wide-band 
approximation E(e) = — iT. 

3.3. Self- consistent Hartree approximation 

The simplest way to close Eq. (|17|) is to perform the truncation by approximating 

((xd,<ft))M{x)((d,rf)). (18) 
Then we obtain immediately for the Green function 

G, "< g > = «"•"'» - rns J w+ ,r ' < 19 » 

Here the quantity (x) remains to be calculated. To this end we compute, respectively, 
the time derivatives of the x-operator, 
3x 

i— = lx,H]=cj p, (20) 



and the p-operator 
.dp 



i— = ujqx + 2Xd^d . (21) 
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Upon combining Eqs. (f2T))) and (J^J) we get 

- ^4 = ufa + 2XLJ dU. (22) 

In the stationary case addressed, Eq. (f2"2")) yields a direct connection between the 
"position" of the vibron and the particle number in the dot: 

(x} = -2— (d f d) =-2—n. (23) 

W ^0 



In view of Eq. p9|) . we finally obtain for the spectral function ([6]) the following self- 
consistent expression , 

2F 

A{e) = - - 2 . (24) 

(e-e~o + 2^n) + P 

This result is equivalent to the one obtained earlier in Rcfs. [24, 25, 20j using alternative 
approaches. The same spectral function can be found if one takes the self-energy in 
Hartree approximation. 

3. 4- Second approximation 

In the first approximation above, the self-consistent Hartree treatment, fluctuations of 
the particle number n and the vibron coordinate x are completely neglected. In order 
to go one step further and estimate possible corrections, we start from the generated 
equations for the second-order Green functions, ((xd, d*)) and ((pd, d^)), 

(e + in - e ) ((xd, d 1 *)) 

= (x) +u;o((pd,d^) + X((x 2 d,d^)) +J2v* k ((xc hk ,S)) , (25) 

i .k 

{e + ir,-i )((pd,d ] )) 

= lj ((xd, dt)) + \((pxd, dt)) + V* k ((pCi >k ,S)) . (26) 

i . k 

The second approximation that we consider here is based on the factorization 

((x 2 d,J)) « (^((ad,^)), 
((p^dt)) » (x)«pd,dt>> + 2«d,dt», 
((xc^ic, d f )) « (x) ((ci,fc, d f )) , 

«pc^,d t ))«0. (27) 

The corresponding set of equations reads 

(e + M7-e -E) ((<*,<**)) = 1 + A((a;d,d t )), 

(e + ^-eo-A^X^dt)) = (x}+Lo ((pd,d^) + (x)Z((d,J)), 

(e + in-eo-X (x)) ((pd, d+» = w ((xd, d f » + 2A ((d, rft)) , (28) 

from which we obtain the second approximation for the Green function, 



[G"-(e)]- x = [G^)]- 1 - ( A (x) + - ) , (29) 



A 2 A n 



A-fl' 
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Figure 2. (Color) Bistable behavior of the level population n in the self-consistent 
Hartree approximation (left) and second approximation (right) as a function of 
gate voltage for different electron- vibron interaction strength A = 3ujq (black), 
A = 4ljo (red), A = huiQ (green), the other parameters are Y = 5u;o, T = 0.25wo, 
r) = and Vb = 0. 



where we introduced A = e + ir\ — e , Q = a-\(x) ' an< ^ * s gi ven by the 
Green function obtained in the Hartree approximation, Eq. (|19j) . After inserting the 
expression (|23|) for the level population (x), Eq. |29|) reads 

[G^e)}- 1 = [G'^e)}- 1 (1 - 2n) • (30) 

We then calculate the spectral function © and the average number of electrons, ©. 
The self-consistent calculation is performed following the chain G r (e) — > A(e) — > n — > 
C(e). 

Before entering into the discussion of the calculated quantities, we consider the 
structure of the Green functions obtained in the two different approximations: 

• In the limit where Q — > the second approximation reduces to the first one. 
Although the EOM method is not a systematic expansion, it tells us that the 
second approximation consistently extends the first one and reproduces it in a 
limiting case. 

• The second term on the right site of Eq. (|3Q[) represents an additional shift with 
respect to Eq. (fT9|) . In the case of very small frequencies ujq, Eq. (|30|) reduces to 
[G r H (e)]~ — (1 — 2n) A 2 a(a-a(.t)) ' involving the first term of a series expansion 
in luq. 

4. Results, discussion and conclusions 

Starting from the expressions derived above within the first (self-consistent Hartree) 
and second approximation we have performed numerical simulations for the average 
population and current in the molecular junction as a function of gate and bias voltage 
for the different parameters A, T and T entering the model Hamiltonian. Furthermore 
we compare the two approximations. The underlying nonlinear equations give rise to 
bistability in the level population thereby enabling memory effects and affecting the 
current. Below we analyze in detail the parameter ranges and conditions for memory 
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Figure 3. (Color) Level population n in the self-consistent Hartree 
approximation (left) and the second approximation (right) as a function of 
normalized bias voltage for different electron-vibron interaction strength A = 3u>o 
(black), A = 4ujo (red), A = 5u>o (green), the other parameters are T = 5ujo, 
T = 0.25oj , V = and V G = 



effects to occur. We focus here on the case of an asymmetric junction, rj = 0, while 
keeping the coupling to the leads symmetric, = Tr = T, to reduce the parameter 
space. 

First we investigate the gate- voltage dependence of the level population. In Fig. [2] 
it is clearly seen that bistability takes place only at larger values of the electron-vibron 
coupling A. The critical values at which bistability occurs and disappears depend 
on the coupling T to the leads and on temperature T; we discuss this parameter 
dependence below. Note that at large values of A the level population of the stable 
states is close to and 1; thus these two memory states are well distinguishable in 
charge. 

A comparison between the left and right panel of Fig. [51 representing the two 
different levels of approximations in the EOM method, shows that the self-consistent 
Hartree treatment underestimates the parameter range where bistability occurs: the 
critical value of A for the occurence of bistable behavior is close to 3wo (see left panel) . 
At this value a bistable regime has already developed for the second approximation. 
This can be partially understood taking into account the additional term appearing 
in the Green function for the second approximation, Eq. (|30j) . This term increases the 
polaron shift thereby enhancing the bistable behaviour in the second approximation. 
Fig. [3] shows the bias-voltage dependence of the level population. In this case 
a qualitative difference arises between electrostatically symmetric (r/ = 0.5) and 
asymmetric (?y = 0) junctions. For asymmetric junctions both states are stable at 
zero voltage, and both charge states are easily accessible. As we showed recently 
|23j . asymmetric junctions are thus favourable, since they exhibit memory effects and 
hysteretic behavior at zero bias, enabling controlled switching upon ramping the bias 
voltage. For symmetric junctions hysteresis is expected only at finite bias voltage 
(noncquilibrium bistability [20]), and hence only a single stable state exists at zero 
bias. Furthermore, at finite voltage the level is only partially occupied by tunneling 
electrons. These two features render the symmetric system less suited for a memory 
setup. 
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Figure 4. (Color) Current normalized to To = Te/4 versus bias voltage for the 
self-consistent Hartree approximation (left) and the second approximation (right) 
for different electron- vibron interaction strength A = 3wo (black), A = 4u>o (red), 
A = 5uq (green), for the same parameters as in Fig. [3] 



We further display in Fig.[4]the current-voltage characteristics, which reflects the 
switching behavior of the system. The characteristic feature is a current jump at the 
bias voltage value where recharging sets in. This behavior can be used to test the 
state of the system and as readout. 

Finally, we depict in Fig. O "phase diagrams" showing the boundaries between 
the parameter regions where bistable memory states (below the boundaries in the 
two panels of Fig. [5|) and single-valued states exist. The left panel of Fig. [5] shows 
the r — T parameter plane. The curve separating single-valued and bistable states 
can be roughly approximated by the condition T + T = c(A), where c(A) is extracted 
to be c(A) m A 1 ' 7 . This means that either thermal or quantum tunneling broadening 
suppresses the hysteresis. The disappearance of the memory effect at high temperature 
is due to enhanced electron tunneling into the higher-energy state. We note that in the 
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Figure 5. (Color) "Phase diagram" depicting the boundaries between parameter 
regimes of bistable memory (below the threshold lines) and single-valued states 
(above threshold). Left panel: F — T parameter plane at different electron- vibron 
interaction strength A = 3loo (black), A = 4u>o (red): ^ = 5^0 (green); right panel: 
A — r paramter plane for T = 0.25iJo- 
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analysis of the temperature dependence the effect of vibronic coordinate fluctuations 
is neglected that can be relevant close to the threshold between single- valued and 
bistable regime. The right panel in Fig. displaying the A 2 — T parameter plane, 
shows on the one hand that bistable behaviour requires, for growing T, increasing 
electron-vibron coupling A. The boundary between the two regimes is approximately 
a straight line in the A 2 — T plane. Hence the condition for finding the memory 
effect is given by T < 0.63A 2 /wo- This clearly shows that for the appearence of the 
memory effect at low temperature, the two energy scales to be compared are the level 
broadening T and the polaron shift A 2 /a>o- A further important conclusion is that the 
memory effect is suppressed for large coupling to the leads. Since, on the other hand, 
larger coupling favours fast information writing and reading and also can additionally 
suppress effects from quantum tunneling between states, the problem arises to find 
optimal parameters for utilizing the memory effect. This will be the subject of future 
work. 

To conclude, we considered a charge-memory effect and switching phenomena 
within a single-level polaron model of a molecular junction in the framework of the 
cquation-of-motion approach to the noncquilibrium Green function theory at different 
levels of approximation. Electrostatically symmetric and asymmetric junctions show 
qualitatively different bistability behavior. In the latter case, controlled switching 
of the molecule is achieved by applying finite voltage pulses. We showed that 
bistability takes place for sufficiently large elelctron-vibron coupling for a wide range 
of further parameters such as molecule-lead coupling and temperature and performed 
a systematic analysis of this parameter dependence by computing phase diagrams for 
the memory effect. We focussed on the investigation of bistability for stationary states, 
leaving the problem of time-dependent fluctuations and fluctuation-induced switching 
as an open question to be investigated in future research. 
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